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FULLY MULTIDIMENSIONAL FLUX-CORRECTED TRANSPORT 


I. INTRODUCTION: FCT DEFINED 


Consider the following system of equations 

w, + f x . = 0 ( 1 ) 

where w and /are vector functions of independent variables xand r. A simple example of such 

a system of equations would be the one dimensional equations of ideal, inviscid fluid flow: 


w = 

p 1 

pv ; 

/ = 

pv 

pv 2 -1- P 


pE) 


p\E + P v 


where p, v, P and E are the fluid density, velocity, hydrostatic pressure and specific total energy 
respectively. 

We shall say that a finite difference approximation to Eq. (1) is in conservation (or "flux") 
form when it can be written in the form 

< + l = <“ & X ,-' |f, + ( i/ 2 ) - ^-(1/2)] (2) 

Here w and / are defined at the spatial grid points x, and temporal grid points t", and Ax, = 


y (x, +1 — x,. The F, +(l/2 ) are called transportive fluxes, and are functions of / at one or 

more of the time levels t". The functional dependence of F on fdefines the integration scheme 
(leapfrog, Lax-Wendroff, Crank Nicholson, donor cell, etc.). 

It is well known that higher order (order 2 and above) schemes for numerically integrat¬ 
ing Eq. (1) suffer from dispersive "ripples" in w, particularly near steep gradients in w. Lower 
order schemes, such as donor cell, Lax-Friedrichs, or high order schemes with a zeroth order 
diffusion added, produce no ripples but suffer from excessive numerical diffusion. Flux- 
corrected transport (FCT) is a technique developed by Boris and Book [1-3] which embodies 
the best of both of the above worlds. In its simplest terms, FCT constructs the net transportive 
flux point by point (non linearly) as a weighted average of a flux computed by a low order 
scheme and a flux computed by a high order scheme. The weighting is done in a manner which 
insures that the high order flux is used to the greatest extent possible without introducing rip- 
’Manuscript submitted May 17, 1978 
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pies (overshoots and undershoots). This weighing procedure is referred to as "flux-correction" 
or "flux-limiting 1 ' for reasons which shall become clear later. The result is a family of transport 
algorithms capable of resolving moving contact discontinuities over 3-4 grid points, and shock 
fronts over 2 grid points, without undershoot or overshoot [1-31. Formally, the procedure is as 
follows: 

1) Compute F,+ (1/2 ), the transportive flux given by some low order scheme guaranteed to 
give monotonic (ripple-free) results for the problem at hand 

2) Compute /\+ (l/2 ), the transportive flux given by some high order scheme 

3) Define the "amidiffusive flux": 

A,+t\m = ^7+0/2) - ^m-U/2) 

4) Compute the updated low order ("transported and diffused") solution: 

*,;</= w »- Ax,’ 1 [f,+{| /2 ) - F[. l (1/2) ] 

5) Limit the /4, +(1/2 ) in a manner such that w ,, + l as computed in step 6 below is free of 
overshoots and undershoots: 

/4 / + (l/2) = C+(l/2> ^-+(l/2)< 0 ^ C+(l/2) ^ 1 

6) Apply the limited antidiffusive fluxes: 

< + ‘ = w ', d - Ax, - ' [>4 i+(1/2) ~ ^-(1/2)] 

The critical step in the above is, of course, step 5 which will be discussed shortly. In the 
absence of the flux limiting step (>f, c +u/ 2 > = -4,+n/ 2 >), w" +] would simply be the time-advanced, 
high order solution. 

We note that this definition of FCT is considerably more general than has been given pre¬ 
viously. 

II. MULTIDIMENSIONAL FLUX-CORRECTED TRANSPORT 

Before proceeding to a discussion of flux limiting, let us see how the procedure given 

above might be implemented in multidimensions. An obvious choice would be to use a 

Strang-type time-splitting procedure [4] when it can be shown that the equations allow such a 
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technique to be used without serious error. Indeed, such a procedure may even be preferable 
from programming and time-step considerations. However, there are many problems for which 
time-splitting produces unacceptable numerical results, among which are those involving 
incompressible or nearly incompressible flow fields. This technique is straightforward and shall 
not be discussed here. Instead we consider as an example the fully two-dimensional set of equa¬ 
tions 


H>, + f X + By “ 0 (3) 

where w, /, and g are vector functions of x, y , and I. In finite difference flux form we have 

= w "j ~ A Vj~j |^'/+(l/2), J - ^—(1/2). j + Gi.j+( 1/2) ~ ^/.2-(l/2)] (4) 

where now w, / and g are defined on spatial grid points x,, yj at time levels i", and A V, j is a 
two dimensional area element centered on grid point (/, j). Now there are two sets of tran¬ 
sport! ve fluxes F and G , and the FCT algorithm proceeds as before: 

1) Compute Ff+ W 2 ).j and G,f- /+( j/ 2) by a low order monotonic scheme 

2) Compute F,+ (1 / 2 ).j and G, H J+(l/2 ) by a high order scheme 

3) Define the antidiffusive fluxes: 

^<+(1/2).; = ^-+(1/2).; ~ ^/+(i/2), / 

^i,J+(\/2) = G, W y + ( 1/2) - y+(l/2) 

4) Compute the low order time advanced solution: 

w !fj = w ?j ~ [^+(i/2), j ~ ^-d/2), j + GuH\m ~ Gu-nn >] 

5) Limit the antidiffusive fluxes 

/i f+(V2).j = -4, + (1/2), j G, + (l/2 ).j 0 ^ G, + ( i/2),y < 1 

^U+(i /2) = -41, y+d/2), y G, y+d/2) 0 < C/.;+(i/ 2 ) ^ 1 

6) Apply the limited antidiffusive fluxes: 

w ’"y l = w kj “ A V<7.; [^z+d/a)./ " Af-(\n),j + 4, C j +( i/2) - <4,fy_d/2)| 

As can be easily seen, implementation of FCT in multidimensions is straightforward with 
the exception of Step 5, an algorithm for which is the subject of this paper. First, however, let 
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us see how flux limiting is presently implemented in one spatial dimension. 


III. FLUX LIMITING IN ONE SPATIAL DIMENSION - 
THE ORIGINAL ALGORITHM 

The original algorithm for flux-limiting in one dimension was given by Boris and Book 
[1]. In our notation it is: 


A,+ (1/2) — -S/-MI/2) max 


0, min \A, + ( 


^+(1/2) w !+2 ~ w !ti\ Ax-, +I , 5 , , +n/2) Iw,"'- tv/fO A.v,_,l 


(5) 


where 


S/ + U/2) - 


+ 1 if -4,+d/2) ^ 0 
— 1 if /4,+d/2) < 0 


The intent of this formula is that antidiffusive fluxes should neither create new extrema, nor 
accentuate already existing extrema, in the transported and diffused solution w"'. That the 


above formula does, in fact, perform precisely this task can be verified by the reader with rela¬ 
tive ease. We shall examine here some of the less obvious properties of this very powerful, yet 
simple, formula. In the process we shall gain insight into which of these properties we shall 
wish to carry over into a multidimensional flux limiter. We first observe that certain quantities 
do not appear in the above formula: 1) - w,"', the first difference of tv"'at the point where 

the antidiffusive flux -4,+ <l/2 ) is evaluated; and 2) antidiffusive fluxes other than A, + (l/2) . This 
last property is the most notable since there are conceivably two fluxes directed into or out of a 
cell. A formula guaranteeing that the two fluxes acting in concert shall not create ripples would 
apparently require a knowledge of both. We shall return to this point momentarily. 

In Figure I we show the eight possible configurations of tv"' in the neighborhood of a 
positive <4,+(i/ 2) ( directed to the right in our diagrams). Configurations 1-4 show the "normal' 1 
situation, with /4,+ (l/2) having the same sign as tv,"', - tv,"' (as might be expected of an 
antidiffusive flux ). We note that if either *v,'+ 2 — tv,'+j or tv,"' — w't \ has a sign opposed to that 
A,+u/ 2 ), as in configurations 2-4, the antidiffusive flux -4, +-r i / 2 > is completely canceled. This, 
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however, is in total agreement with the stated intent of Eq. (5) since otherwise configuration 2 
would allow accentuation of an existing maximum, configuration 3 accentuation of an existing 
minimum, and configuration 4 accentuation of both. In the remaining configuration 1, the flux 
limiter (5) will reduce the magnitude of A,+{ U2 ) sufficie 1 '"-' to guarantee that neither a max¬ 
imum at grid point / + 1 nor a minimum at grid point / will be formed, again in precise agree¬ 
ment with its stated intent. 

Configurations 5-8 are identical to configurations 1-4, respectively, except that sign of m>,+, 
— w,"'has been reversed (The "antidiffusive fluxes" are now directed down the gradient in w' J ). 
Since the sign of - w'f does not enter into the flux correction formula (5), the results of 
the formula are identical to those for the previous four cases: the antidiffusive fluxes are can¬ 
celed for configurations 6-8 and limited in configuration 5 to the extent necessary to prevent a 
new maximum at grid point / + 1 and a new minimum at grid point /. Examination of 
configurations 6-8 reveals that /f, + (1/2 ) actually presented no hazard insofar as extrema creation 
or enhancement (at least in moderation). Certainly there was no cause for completely cancel¬ 
ing the flux. Even in configuration 5 the flux may have been limited to a greater extent than 
necessary. At first it would seem that configurations 5-8 represent errors introduced by the 
simplicity of the flux limiting formula (5). However, extensive tests by this author indicate 
that in the relatively rare instances in which configurations 5-8 occur in practice, the "errors" 
introduced by Eq. (5) represent, in fact, the correct action to take in terms of producing accu¬ 
rate profiles in w" +l . More importantly, they represent the mechanism by which Lq. (5) can 
guarantee that ripples are not formed under any circumstances, as we shall see presently. 

Consider two antidiffusive fluxes, acting in concert, attempting to produce or accentuate 
an extremum. We therefore have /f, +n/2) and -4,-n/j) either both directed toward, or both 
directed away from grid point i. We see from Figure 1 that, in general, an antidiffusive flux 
directed opposite to the gradient in w ,rf will be completely canceled. Therefore the only cases of 
fluxes acting in concert that we need be concerned with are those where two adjacent fluxes are 
both parallel to the local gradients in w' d . These are precisely the cases of already existing 

, si 

\j */ 
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extrema, in which case both fluxes will be canceled (as in configurations 2-4). This is the reason 
that Eq. (5) needs no information on any antidiffusive flux other than A l+(i/2) . 

In Figure 2 we see that the above-mentioned assumptions regarding antidiffusive fluxes 
acting in concert break down completely in multidimensions. It is possible in more than one 
dimension for more than one antidiffusive flux to be directed into or out of a cell, all of these 
fluxes being directed parallel to the local gradient in w " 1 , without that cell being an already 
existing extremum. Therefore the problem of dealing with multiple antidiffusive fluxes acting in 
concert cannot be avoided by simply canceling all fluxes antiparallel to the local gradient in w' d . 
It is clear then that any formula which purports to perform flux limiting in more than one 
dimension without resort to time splitting must contain information about antidiffusive fluxes 
other than the one being limited. 

IV. FLUX LIMITING IN ONE SPATIAL DIMENSION - 
AN ALTERNATIVE ALGORITHM 

We describe here in one spatial dimension an alternative flux limiting algorithm which 
generalizes easily to multidimensions and which, even in one dimension, exhibits a superiority 
to the limiter described in the previous section (Eq, (5)) with regard to peaked profiles. 

Referring to Figure 3, we seek to limit the antidiffusive flux A, +{xm such that 

^M-d/2) = C + tl/2) 'L + tl/2)' 0 ^ C + O/2) ^ 1 (6) 

and such that /f,+o/ 2 ) acting in concert with A^ (l/2 ) will not allow 

w"*' = W,' d - A.X, 1 J/f, C + (j/2) - -4, C (1/2)] 

to exceed some maximum value w, max nor fall below some minimum value tv, m ' n . We leave the 
determination of w, max and w, mm until later. 

We define three quantities: 

P, + = the sum of all antidiffusive fluxes into grid point / 


= max (0, -4 ,-<i/ 2>) ~ min (0, -4, + u/2>) 

(7) 

Q * = ( H', max — w'^) \x, 

(8) 

I mind. Q, + /P, + ) if P, + > 0 

R ‘ “ 1 o if p; = 0 

(9) 
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Assuming that w, miX > w' d (it must be), all three of the above quantities are positive and R, + 
represents the least upper bound on the fraction which must multiply all antidiffusive fluxes into 
grid point / to guarantee no overshoot at grid point i. 

Similarly we define three corresponding quantities: 

p- = the sum of all antidiffusive fluxes away from grid point i 

= max (0, /t, +u/2 )) - min (0, A ,- U/2 >) (10) 

Q-- (wf- w, m ' n ) bx, (11) 

|min(l, Qr/Pf ) >f P, > ol ( 12 ) 

R > ~ | 0 if Pf = 0 J 

Again assuming that w, mi " < w* we find that R~ represents the least upper bound on the 
fraction which must multiply all antidiffusive fluxes away from grid point / to guarantee no un¬ 
dershoot at grid point /. 

Finally we observe that all antidiffusive fluxes are directed away from one grid point and 
into an adjacent one. Limiting will therefore take place with respect to undershoots for the 
former and with respect to overshoots for the latter. A guarantee that neither event comes to 
pass demands our taking a minimum: 

min (/?,-+ 1 , R,~) if ^i+a/ 2 ) ^ 0 
C, + (l/2) - m j n (R+ l /?,+]) if A ,+(i/2) < 0 

Furthermore, we shall call upon our previously described experience with the original flux 
limiter and set 

A, + o/2) = 0 if A, +u/2 ) (wAi - <0 (14) 

and either A i+ o/ 2)(* v /+2 ~ w !+0 < ® 

or A,^ m) (w! d - W&) <0 

In practice the effect of Eq. (14) is minimal and is primarily cosmetic in nature. This is 
because cases of antidiffusive fluxes directed down gradients in w’ d are rare, and even when 
they occur usually involve flux magnitudes that are small compared to adjacent fluxes. If Eq. (14) 
is used, it should be applied before Eq. (6) through (13). 
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We come now to a determination of the quantities w’, max and w, mm in Eqs. (8) and (11). 
A safe choice is 

w max = max w! d , w u {,) (15) 

n min = min ( w' d |, w' J , ,) (16) 

This choice will produce results identical with those of Eq. (5) in one dimension, including the 

occurrence of the "clipping" phenomenon to be mentioned shortly. 

A better choice is: 


w, = max ( w", w ) 


; = max (w,"_ |, w a , wj +,) 

(17) 

w h = min (w”, w," 1 ) 


1 = min (wi, R’*, w/Vj) 

(18) 


This choice allows us to look back to the previous time step for upper and lower bounds on 

w," +1 . 

It is clear that these two methods of determining w, max and w, min represent only a small 
subset of possible methods. The alternative flux limiter described in equations (6) through (14) 
admits of any physically motivated upper and lower bound on w," 1 supplied by the user, intro¬ 
ducing a flexibility unavailable with the original flux limiter (5). However, with the exception 
of one example in the next section (which shows graphically the potential power of this flexibil¬ 
ity), we shall henceforth use Eq. (17) and (18) to evaluate w min and tv max in one dimension. 

V. COMPUTATIONAL EXAMPLES - ONE DIMENSION 

We consider one dimensional passive convection in a constant velocity field. We have Eq. 
(1) with w = f) and / — pv with v = constant. We choose our transport algorithm to be that 
given in [3] for LPE Shasta. On the standard square wave tests we find that our results for the 
original flux limiter (5) and for the alternative flux limiter (6) through (14) are identical to 
within round-off (the same is true for traveling shock waves in the coupled one dimensional 
equations of ideal inviscid fluid flow). To find differences between the limiters in one dimen¬ 
sion we must look to passive convection of peaked profiles. We choose the problem given by 

Forester [5], a guassian of half-width 2A.v. In Figure 4 we show the results after 600 iterations 
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for the trivial case v = 0. On the left we see the fimilar "clipping" phenomenon with the origi¬ 
nal flux limiter, caused by a zeroth order diffusion term in the low order portion of the LPE 
Shasta algorithm. This diffusion term causes the peak in w'^to be smaller than the peak in w,", 
leaving the original flux limiter (5) with no way of resurrecting the original peak. This process 
occurs repeatedly, eventually leaving the characteristic three point top. The alternative flux 
limiter, shown on the right, "remembers" the old value of the peak and is able to resurrect it 
each time step. 

In Figure 5 we show the same problem after 600 iterations but this time for e = v At/ Ax 
= 0.1. We see that clipping occurs with both flux limiters, but to a lesser extent with the alter¬ 
native flux limiter (6) through (14). 

At this point we removed the flux limiter entirely and again ran the problem 600 itera¬ 
tions with e = 0.1. The results convinced us that the amplitude and phase properties of the 
high order portion of LPE Shasta were incapable of resolving the high wave numbers of which 
the gaussian is composed. Consequently it was decided to switch to a higher order algorithm, a 
leapfrog-trapezoidal transport algorithm which uses eighth order spatial differences. The algo¬ 
rithm is, then, second order accurate in time and eighth order in space, with an amplification 
factor that is effectively unity across the entire Fourier spectrum, and phase properties consider¬ 
ably better than those of fourth order algorithms. The leapfrog portion of this algorithm is 
identical to the eighth order Kreiss-Oliger scheme [6]. A fourth order version of this same 
algorithm was used later in the two-dimensional solid body rotation tests. We ran the gaussian 
test problem again 600 iterations with « = 0.1 with no flux limiter and were convinced that the 
algorithm did indeed have the resolving power necessary to do the problem. A low order 
scheme, donor cell plus a zeroth order diffusion term with coefficient -j-, was added to com- 

O 

plete the FCT algorithm, which we dub 2-8 leapfrog-trapezoidal. A more detailed description 
of this algorithm is found in the appendix. 

Figure 6 shows the results of 2-8 leapfrog trapezoidal run 600 iterations with e = 0.1 with 
both the original and alternative flux limiters. The results are better than those in Figure 5, 
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and again the alternative flux limiter proves superior, but nontheless disappointing. The clip¬ 
ping would appear to be due entirely to the flux limiters, not to the phase or amplitude proper¬ 
ties of the high order scheme. 

A careful examination of exactly what happens to a one point peak in a finite difference 
code reveals the real source of the above problem. Consider a profile with a local peak at grid 
point i in passive convection at constant velocity > 0. At each succeeding time step the func¬ 
tion value at grid point / will decrease and that at / + 1 will increase (Figure 7). Eventually 
they will both reach some intermediate value, and the actual original peak value will not appear 
anywhere on the grid, since it’s position now lies midway between two grid points. At this 
point even the new flux limiter (6) through (14), (17) and (18), has lost the information it 
needs to allow the peak to be resurrected in suceeding time steps, and will "clip" the new peak 
at grid point i + 1 as it tries to form, based on the assumpi : on that it is, in fact, an overshoot. 
The effect is magnified, since the clipping itself introduces phase errors in succeeding iterations, 
the net result being the profiles dipicted in Figure 6. 

It is clear, then, that if we are to successfully treat a one-point extremum within the con¬ 
text of FCT we must use information other than just the grid point values themselves. In what 
follows we shall utilize the flexibility of the alternative flux limiter to use as tv, max and tv," 11 " any 
values that we choose. In Figure 8 we show a possible way of extracting information about 
extrema which do not lie exactly on a grid point at the time. Basically we define tvfflw) to be 
the tv value at the intersection of the line segments formed by connecting the point (x,_i, w,_,) 
with (x„ tv/*) and the point (x, +1 , w/*,) with (x, +2 , tv/? 2 ). If the x coordinate of this intersec¬ 
tion lies between x, and x, + 1 , then we consider this tv, p + e f, k /2) to be an allowable tv max or tv m,n for 
either tv," + 1 or tv/VV- We now have 

w, a = max (tv,", tv/*) 

tv, max = max (tv/L,, tv, D , w,° + ,, w,^f, k / 2)> *v ( ^5f, k /2) ) (19) 

tv* = min (tv,", tv/*) 

tv, min = min (tv*_,, tv* tv,* +1> tv ( ?tf /2) . w, p ', a l k /2) ) (20) 

Equations (19) and (20) together with equations (6) through (14) now determine the 
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alternative flux limiter (for this section only). 

In Figure 9 we show the results of using Equations (19) and (20) to determine w, max and 
w, mm on the gaussian test problem run 600 iterations with e = 0.1. Clearly the problem has 
been solved — we recover the gaussian profile with no dispersive ripples and minimal loss in 
amplitude. We have not performed this test merely to show the power of the extrapolation 
technique just described to determine w, max and w, mm . Rather this calculation serves to show the 
power of using information other than that available on the one dimensionhI grid. In multidimen¬ 
sional flux limiting, this information comes from the other coordinate directions, as we shall 
see. 

VI. FLUX LIMITING IN MULTIDIMENSIONS 

The alternative flux limiting algorithm presented in section IV generalizes trivially to any 
number of dimensions. For the sake of completeness we present here the algorithm for two 
spatial dimensions. 

Referring to Figure 10, we seek to limit the antidiffusive fluxes A l+(U2 ),j and /4, iJ+( i/ 2 , 
such that 

A r+i\m.j = C/-HI/ 2 ).j ‘ 4 - + (l/ 2 ).j 0 ^ C /+(l/ 2),2 ^ 1 
Afj+U/l) = C.j-Ml/2) A, y+(l/2) 0 < C, y+(l/2) < 1 

and such that AfL {m) j , Af J+{ 1/2 ), and ^, c ,_n/ 2 > acting in concert shall not cause 

w n+\ = w rd f _ ^ y-j |/|C (1/2) J - Af_ ivi ),j + -4, C >+( 1/2) “ ^- C j-(l/2)] 

to exceed some maximum value w, max nor fall below some minimum value tv™ 11 . 

Again we compute six quantities completely analogous to those computed in Eq. (7) 
through (12): 

P+j = the sum of all antidiffusive fluxes into grid point (/j j) (7 ) 

- max (0, A,-( 1/2 ),j) - min "0, /4,+n/2),y) 

+ max (0, /<,.y-o/ 2 )) - min '0, /(, j +( i /2 )) 

Q+j _ (h,™>* _ w «j) A V, j (8') 

(mind, Q+j/Pfj if Pfj > 0 
Rj o if/ft-0 
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Pjj — the sum of all antidiffusive fluxes away from grid point O', J) 


- max (0, -4/+0/2). j) - min (0, 

+ max (0, A IJ+(m) ) - min (0, y-<i/ 2 >) 

QCj “ ( w u - < 7 ) AV U 

(mind, Q-j/P-j) if P~j > 0| 

"“I 0 if Kj “ oj 

Equation (13) becomes 


C/+(l/2), j 
C i. J+il/2) 

while Eq. (14) becomes 


jmin (RjX\,j, R,j ) if A i+ ( X / 2 ) j ^ 0 
|min (R+j, Rj+i.j) if A,+ xX / 2 )j < 0 

jmin (R*j+i, Rfj) if <4,.y+<i/j) ^ 0 
Jmin (R, + j, R~j+\ ) if A liJ+ ( |/ 2 ) < 0 


( 10 ') 

on 

( 12 ') 


(13') 


^/+<i/2). y “ 0 if ^,+(1/2), > (W/+I .J - »!?j) < o 

and either A, HU2)iJ (w% 2 j - w#,,) < 0 
or A l+(m) ,j {w'Aj - w\i X J ) < 0 
Ai.j+i i/ 2 ) “ 0 if A iiJ+i i/ 2 ) (tv'< +1 - wj d j) < 0 
and either A iJ+lU2) (,w'?j+ 2 ~ w u+i ) < 0 



or ^/.y+O/2) 

(w' d j - w! d j_ x ) < 0 

(14') 

and Eq. (17) and (18) become 





w tj — 

max (w'j, wfj) 


yy .rn®x — 

max (w,°_i j 

■ w u> w f+\.i- <j- 1. <y+i> 

(17') 


wj’j - 

min(w,^, wj d j) 

(18') 

min s 

min(w, A _ 1 j, 

"Up **'/+!. 7, H *V A y +1 ) 


Again, the effect of Eq. (14') is minimal, but if it is used it should be applied before Eq. 


(6') through (13'). Note that our search for wff x and w, m j n now extends over both coordinate 
directions. Where finite gradients exist in both directions, this procedure will allow us to stop 
the clipping phenomenon in regions where a peak exists with respect to one coordinate direc¬ 
tion but not in the other, as we shall see in the next section. 


VII. COMPUTATIONAL EXAMPLES - TWO DIMENSIONS 

We choose as our two dimensional test problem that of solid body rotation. That is, we 
have Eq. (3) with / - w y x , g - m,, v, - ~(l(y - y a ), and v, - (l(x - x„). Here fl is the 
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(constant) angular velocity in radians/sec and ( x„, y„) is the axis of rotation. The configuration 
is shown in Figure 11. The computational grid is 100 x 100 cells, Ax = Ay, with counterclock¬ 
wise rotation taking place about grid point (50, 50). Centered at grid point (50, 75) is a cylinder 
of radius 15 grid points, through which a slot has been cut of width 5 grid points. The time 
step and rotational speed are chosen such that 628 time steps will effect one complete revolu¬ 
tion of the cylinder about the central point. A perspective view of the initial conditions is 
shown in Figure 12. 

Our high order scheme for the following tests is a fully two dimensional, fourth order in 
space, second order in time leapfrog-trapezoidal scheme, the leapfrog step of which is a two 
dimensional fourth order Kieiss-Oliger scheme [6]. The low order scheme is simply two 
dimensional donor cell plus a two dimensional zeroth order diffusion term with diffusion 

coefficient -j. A more detailed description of this algorithm is found in the appendix. 


We wish to emphasize that the only difference between calculations in the following com¬ 
parisons is in the flux limiting stage itself. The high order fluxes, low order fluxes, and hence 
the (unlimited) antidiffusive fluxes are all computed in the full two dimensions without using 
time-splitting. In each case we are comparing the fully two dimensional flux limiter given by 

Eq (6') through (14'), (17') and (18 ) with a time split application of Eq. (5). Time splitting is 

the only way that Eq. (5) may be utilized in a multidimensional problem. Note that in the 
latter case we are not time splitting the entire transport operator, but only the flux limiter (5). 
In this way we are testing only the limiters themselves. 

In Figure 13 we show a perspective view of the two calculations after revolution (157 

iterations). Figure 14 presents a comparison of the results of the two calculations for one full 
revolution (628 cycles). Two features are obvious. The first is a much greater filling-in of the 
slot with the time split Eq. (5) than with the fully two dimensional flux limiter. The second is 

the loss of the bridge connecting the two halves of the cylinder in the case of the time-split 

application ot Eq. (5). Less obvious is the lack of dipping of the peaked profiles defining the 






front surface of the cylinder for the case of the fully multidimensional limiter. Clearly this is 
due to the fact that the multidimensional flux limiter can look in both directions to determine 
whether or not a genuine maximum exists. Note that there are two factors working in favor of 
the fully multidimensional flux limiter: 1) the ability to look in both directions to find minima 
and maxima, as just mentioned; and 2) the ability to scan both w" y and wfj to find maxima and 
minima. Both of these factors are responsible for the improved profiles. 

VIII. THE STRIATIONS CODE - A TWO DIMENSIONAL INCOMPRESSIBLE 
FLUID CODE USING FULLY MULTIDIMENSIONAL FCT 

A two dimensional (x - y) plasma cloud initialized in a region of constant magnetic field 
B„ directed along the z axis, with an externally imposed electric field E 0 directed along the x 
axis will tend to drift in the E 0 x B„ direction (along the negative y axis) (see Figure 15). If 
the ion-neutral collision frequency is finite, Pedersen conductivity effects will produce polariza¬ 
tion fields which tend to shield the inner (more dense) regions of the cloud from E„, causing 
this inner portion of the cloud to drift more slowly than the outer portions of the cloud. This 
results in a steepening of gradients on the back side of the cloud. Arguments similar to those 
above, applied to infinitesimal perturbations imposed upon this back side gradient, show that 
the back side of the cloud is physically unstable to perturbations along x. For a detailed 
description of this problem, see [7], 

The equations of motion for the electron fluid are: 

(dNJdt) +V ± (N e \ e ) =0 (21) 

V A - (A^V) - E„ ■ V ± N t (22) 

V P = — (c/B 0 ) VjV x £ (23) 

Here N e , V„ and 'P are the electron density, electron velocity and perturbation electric 

d d 

field potential respectively, and V ± is the two dimensional divergence operator x — + y — . 

The magnitudes of B„ and E 0 are 0.5 gauss and 5 millivolts per meter respectively. Our rest 
frame here is that of the (c/B„) E„ x £ velocity. A few trivial vector identities will convince the 
reader that V^V, - 0, meaning that the electrons move incompressibly. Clearly time-splitting 
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the transport operator would be disastrous here, and a fully two dimensional scheme is re¬ 
quired. 

Eq. (22) is solved for ^ using an elliptic solver, and Eq. (23) then yields the electron 
velocity field. We then utilize exactly the same multidimensional FCT transport algorithm used 
in the previous section for solid body rotation to integrate Eq. (21) in time. 

Our computational mesh consists of 40 grid points in the x direction, 160 grid points in 
the y direction, periodic boundary conditions in both directions, and \x = A y — 0.31 km. Our 
"cloud" consists of a 1-D gaussian: 

N e (x,y) = N 0 (1 + lOe (> >0> /64 ) 

where N 0 is the ambient background electron density and y 0 is the spatial center of the gaussian 
distribution. Superimposed upon this distribution is a random x-dependent perturbation with a 
maximum amplitude of 3 percent. 

Figures 16-20 show isodensity contours of N e /N 0 for the above configuration at various 
times in the integration. It is seen that, as expected, the back of the cloud (the upper half in 
the plots) is unstable, growing linearly in the very early stages of development. Non-linear 
effects soon enter the physics, however, as each striation successively bifurcates, producing 
smaller and smaller scale structures, in agreement with the results of the ionospheric barium 
cloud releases which we are attempting to model. Two points which bear on the numerics 
should be noted: 1) the intense gradients dictated by the physics are not diffused away, nor do 
there appear in the problem any of the "ripples" associated with numerical dispersion which 
normally appear when steep gradients try to form; 2) precisely because we did not have to 
resort to time-splitting, none of the usual time splitting phenomena, such as temporal density 
oscillations ans spurious density values, are evident. 

CONCLUSIONS 

We have shown that the algorithm presented in Eq. (6') through (14'), (17 ) and (18 ) 
does, in fact, represent a workable multidimensional flux limiter. In addition, due to the flexi¬ 
bility in determining overshoot and undershoot criteria inherent in the method, the algorithm 
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produces results which are consistently equai or superior to those produced using a lime-spli; 
version of the original flux limiter (5), at least for the admittedly limited class of problems 
presented here. 

For multidimensional problems where time splitting is unacceptable, or for problems 
where the "clipping” phenomenon associated with the original flux limiter (5) is a serious prob¬ 
lem, the new algorithm presented here represents the only way that FCT may be implemented. 
For these problems the choice is clear, for there is only one option. Yet even in situations 
where the constraints mentioned above do not apply, benefit may be gained by implementing 
the new algorithm rather than time splitting Eq. (5). We do not yet have enough experience to 
give any guidelines, and can only ask the prospective user to try the method. 

Certainly the possiblities for modifying the basic scheme are endless. One could, for 
instance, limit the antidiffusive fluxes only with respect to maxima, or to minima; or he could 
limit the fluxes sequentially for maxima and minima, rather than limiting maxima and minima 
simultaneously in the manner presented here (this last procedure will introduce an asymmetry 
between the treatment of maxima and minima which may or may not be desirable). Even 
within time-split codes there are possibilities. One could time split the one dimensional form of 
the new algorithm rather than time splitting Eq. (5); or fully multidimensional flux limiting 
could be performed at the end of each sweep of a time-split scheme. 

On NRL’s Texas Instruments ASC computer, the calculations presented in section VII 
required 93 seconds and 125 seconds of CPU time for the time-split and fully multidimensional 
cases respectively, a cost penalty of slightly more than 30% for the multidimensional limiter. 
Of course this extra cost is highly problem dependent. For instance the striations code 
described in section VIII spends 80% of its time solving Eq. (22), making the net cost penalty 
of fully multidimensional flux limiting only a few percent. 
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Sixth order : 


+ (1/2) ~ |J(/, + i + /,) - ^r(/, +2 + /, i) 
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+ ^(/, +J + /,- 2 ) 


Eighth order : 


+ -no </•-+/■ .> 

The above fourth and eighth order forms are used as the high order (luxes in the main 
body of this paper. 

The low order flux of the leapfrog-trapezoidal FCT schemes is simply donor cell plus 
a zeroth order diffusive flux with coefficient The donor cell algorithm requires that 
/ = vw, where v is a convective velocity. Specifically, 


^ + (1/2) “ v / + (l/2) w '/+ f (t/2) g ( x / + i x <) ( M '/+t 


where 


+ (1/2) = T (v - + V ' + l) 


w" if 


^ 0 


W ' +(V2) -|< +I if v l+(1/2) < 0 

n w"~ ] for leapfrog step 
w u = 

' for trapezoidal step 

A detailed description and analysis of these and other high order FCT algorithms will 
be discussed in a forthcoming report by the present author. 
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—•DENOTES DIRECTION OF A i + y , 



i-1 i i+1 i+2 i-1 i i+1 i+2 i-1 i i+1 i+2 i-1 i i+1i+2 


x 

Fig. 1. The eight possible configurations of the transported and 
diffused solution w ,d in the neighborhood of a positive (rightward- 
directed) antidiffusive flux A, +t i/ 2) . Note that configurations 1 through 
4 differ from configurations 5 through 8 only in the sign of the quantity 
(*•#, - <*>. 



i-1 i i+1 

x 


Fig. 2. Perspective view of a two dimensional profile of the tran¬ 
sported and diffused solution w"\ showing the four possible 
anitdiffusive fluxes affecting the grid point (/, j), the directions of 
which are indicated by arrows. Note that all of the fluxes are paral¬ 
lel to the local gradient in w' d (as "antidiffusive" fluxes might be 
expected to be), and that is not an extremum. This situation is 
impossible in one dimension, and it is precisely this impossibility 
which allows fluxes to be limited without regard to neighboring 
fluxes (see text). In two or more spatial dimensions a flux-limiting 
formula must take into account effects due io multiple fluxes acting 
in concert. 
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i-V* 


w 


td 


• A i+J4 


^ _ 


—I—_I_I- 

i-1 i i + 1 


x 

Fig. 3. One dimensional profile of the transported and diffused 
profile w' d , showing the two antidiffusive fluxes /!,+,,/ 2 > and 
A,_ d/ 2 ) whose collective effect must be taken into account with 
respect to overshoots and undershoots in the final value of 



i 


Fig. 4. Comparison of old and new flux limiters on narrow gaussian 
profile in passive convection for the trivial case of a vanishing velo¬ 
city field. The transport algorithm is LPE SHASTA. Note the 
"clipping" phenomenon associated with the old limiter. 
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Fig. 5. Same comparison as in Fig. 4 except that the 
velocity fluid is now finite. The profile has been con- 
vected through 60 grid points. Note the reduced clip¬ 
ping with the new flux limiter. 



Fig. 6. Same comparison as in Fig. 5, but with a more accurate tran¬ 
sport algorithm (2-8 leapfrog-trapezoidal). Again note the reduced clip¬ 
ping with the new flux limiter. 
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X 


FiR. 7. Time sequence of profiles produced by a "perfect convection 
scheme acting on the variable p with e = 0.2. The actual analytic 
p'ofile is shown as a solid line, and the grid point values are shown as 
dots. Note that at time t 0 + 4At a grid point value at (/ + 1) has 
been generated which is higher than any grid point value at the previ¬ 
ous time step. This is the reason that even the new flux limiter, 
using Eq. (17) and (18) for vv max and w mm , must still "clip". 



Fig. 8. A possible scheme for extracting information about extrema 
which exist between grid points at a given point in time. An extremum 
is assumed to exist between grid points / and i + 1 if the intersection of 
the right and left sided extrapolations of w" has an x coordinate 
between x, and x, +l . The w coordinate of the intersection is then used 
in the computation of w max and w mm (see text). 
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Fig. 9. Same as Fig. 6, except that Eq. (19) and (20), which utilize the 
w peak computation illustrated in Fig. 8., are used to compute w max and 
w m,n in the new flux limiter. Values for the old flux limiter, since they 
are identical to those shown in Fig. 6, are not shown. Note that the 
clipping has been virtually eliminated. 


w td PROFILE: 

j+1 I- 


Jf 


i-i 


i-1 


M-%, j 


A i-»- Va.j 

A i.j-y 2 


i+1 


X 

Fig. 10. Two dimensional profile of the transported and 
diffused values w"', showing the four antidiffusive 
fluxes /f,+(i/ 2 )./, A, -d/ 2 ). /. 1 / 2 ). and A, /- 0 / 2 ) 

whose collective effect must be taken into account with 
respect to overshoots and undershoots in the final value 
of A perspective view of a similar profile is 

shown in Fig. 2. 
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Fig. 11. Schematic representation of two dimensional solid 
body rotation problem. Initially w inside the cut-out 
cylinder is 3.0, while outside w = 1.0. The rotational speed 
is such that one full revolution is effected in 628 cycles. 
The width of the gap separating the two halves of the 
cylinder, as well as the maximum extent of the "bridge" 
connecting the two halves, is 5 cells. 
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Fig. 12. Perspective view of initial conditions for 
the two dimensional solid body rotation problem. 
Note that only a 50 x 50 portion of the mesh cen¬ 
tered on the cylinder is displayed. Grid points 
inside the cylinder have w , = 3.0. All others 
have h\ , = 1.0. 
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Fig. 13. Comparison of perspective views 
of the w profile after 157 iterations (1/4 
revolution) with both the old and new flux 
limiters. The perspective view has been 
rotated with the cylinder, so that direct 
comparison with Fig. 12 can be made. 
Again we plot only the 50 x 50 grid cen¬ 
tered on the analytic center of the cylinder. 
Features to compare are the filling-in of the 
gap, erosion of the "bridge", and the rela¬ 
tive sharpness of the profiles defining the 
front surface of the cylinder. 


NEW LIMITER 



Fig. 14. Same as Fig. 14, but after 628 
iterations (one full revolution). Again note 
decreased diffusion with new flux limiter. 
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PLASMA DENSITY 
ENHANCEMENT 



Fig. 15. Schematic representation of the development 
of a plasma cloud (plasma density increasing toward the 
center) in crossed electric and magnetic fields. Super¬ 
imposed on the bulk E„ x B„ motion is a steepening of 
the rearward side of the cloud. This same side is physi¬ 
cally unstable to small perturbations. 
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Fig. 16. Isodensity contours of 
plasma density at t = 0 sec. The ini¬ 
tial destribution for NJN„ is a gaus- 
sian in y , centered at y — 12.1 km, 
plus a small random perturbation in x. 
Contours are drawn for N,./N„ = 1.5, 
3.5, 5.5, 7.5 and 9.5. The area 
between every other contour line is 
cross-hatched. Only 120 of the 160 
cells actually used in the y direction 
are displayed. Boundary conditions 
are periodic in both directions. In our 
plot B„ is toward the reader, and E„ is 
directed toward the right, and we 
have placed outselves in a frame 
moving with the (c/|fl„| 2 ) E„x B„ 
velocity. The upper portion of the 
gaussian is physically unstable to per¬ 
turbations, while the lower half is 
(linearly) stable. 
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T = 3 0 4 SEC 



Fig. 19. Same as Fig. 16, but for / = 
304 sec. Development is fully non¬ 
linear, as the intense gradients and 
associated high Fourier wave numbers 
become apparent. 


T = 4 0 7 SEC 



Fig. 20. Same as Fig. 16, but for t = 
407 sec. Several plasma bifurcations 
are apparent, in agreement with the 
experimental results from ionospheric 
barium cloud releases, and we have 
maximum to minimum density varia¬ 
tions resolved over only 2 cells. 
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